###########################################################
### Disorganized Political Violence:                    ###
### A Demonstration Case of Temperature and Insurgency  ###
### Replication Code:                                   ###
### Temperature Deviation Results:                      ### 
### Country-Wide Robustness Checks                      ###
### Andrew Shaver, Alex Bollfrass                       ###
### Date: 07/04/2022                                    ###
###########################################################

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load required packages:
library(broom)
library(data.table)
library(radiant.model)
library(lfe)
library(ggplot2)
library(ggpubr)
library(stargazer)
library(raster)
library(rgeos)
library(insol)

packageVersion("broom")
# ‘0.7.12’
packageVersion("data.table")
# ‘1.14.2’
packageVersion("radiant.model")
# ‘1.4.4’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("ggplot2")
# ‘3.3.6’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("stargazer")
# ‘5.2.3’
packageVersion("raster")
# ‘3.5.15’
packageVersion("rgeos")
# ‘0.5.9’
packageVersion("insol")
# ‘1.2.2’

# Load datasets:
panel_af <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/af_panel.Rds")
panel_iq <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iq_panel.Rds")

# Afghanistan boundaries shapefile:
# These .shp files can be accessed at: https://esoc.princeton.edu/data/administrative-boundaries-398-districts
afghanistan_districts <- shapefile("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Raw Data/AFG_district_398/district398.shp")
centroid <- gCentroid(afghanistan_districts, byid = T)

date_seq <- as.data.frame(unique(panel_af$date))
colnames(date_seq) = "date"
date_seq$year <- year(date_seq$date)
date_seq$month <- month(date_seq$date)
date_seq$day <- day(date_seq$date)

temporary <- list()
for(i in 1:length(centroid)){
  temp_matrix <- matrix(NA, length(unique(panel_af$date)), 3)
  for(j in 1:length(unique(panel_af$date))){
    temp_matrix[j, ] <- daylength(centroid$x[i], centroid$y[i], JDymd(date_seq$year[j], date_seq$month[j], date_seq$day[j]), 3)
  }
  # Now, convert to a dataframe and add the district name:
  temp_df <- as.data.frame(temp_matrix)
  temp_df$district <- afghanistan_districts@data$DISTID[i]
  # Then, create date, week and month variables:
  temp_df$date <- unique(panel_af$date)
  temp_df$week <- floor_date(temp_df$date, "weeks")
  temp_df$month <- floor_date(temp_df$date, "months")
  temporary[[i]] <- temp_df
}
daylength <- rbindlist(temporary)

colnames(daylength)[3:4] <- c("daylength", "DISTID")
daylength$DISTID <- as.character(daylength$DISTID)

panel_af <- plyr::join(panel_af, daylength[, c("daylength", "DISTID", "date")])

# Iraq boundaries shapefile:
# These .shp files can be accessed at: https://esoc.princeton.edu/data/administrative-boundaries-districts
iraq_districts <- shapefile("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Raw Data/Iraq Districts/iraq_districts.shp")
centroid <- gCentroid(iraq_districts, byid = T)

date_seq <- as.data.frame(unique(panel_iq$date))
colnames(date_seq) = "date"
date_seq$year <- year(date_seq$date)
date_seq$month <- month(date_seq$date)
date_seq$day <- day(date_seq$date)

temporary <- list()
for(i in 1:length(centroid)){
  temp_matrix <- matrix(NA, length(unique(panel_iq$date)), 3)
  for(j in 1:length(unique(panel_iq$date))){
    temp_matrix[j, ] <- daylength(centroid$x[i], centroid$y[i], JDymd(date_seq$year[j], date_seq$month[j], date_seq$day[j]), 3)
  }
  # Now, convert to a dataframe and add the district name:
  temp_df <- as.data.frame(temp_matrix)
  temp_df$district <- iraq_districts@data$ADM3NAME[i]
  # Then, create date, week and month variables:
  temp_df$date <- unique(panel_iq$date)
  temp_df$week <- floor_date(temp_df$date, "weeks")
  temp_df$month <- floor_date(temp_df$date, "months")
  temporary[[i]] <- temp_df
}
daylength <- rbindlist(temporary)

colnames(daylength)[3:4] <- c("daylength", "ADM3NAME")
panel_iq$ADM3NAME <- as.character(panel_iq$ADM3NAME)

panel_iq <- plyr::join(panel_iq, daylength[, c("daylength", "ADM3NAME", "date")])

# First, generate descriptive statistics:

# PRODUCES TABLES 4 AND 5:
stargazer(panel_af[, c("temperature", "precipitation", "df", "idf", "ied")])
stargazer(panel_iq[, c("temperature", "precipitation", "unconstrained1", "constrained1", "ied")])

# Next, move to temperature deviation analysis:

# Begin by setting range of temperatures:
temperature_range <- 75:105

# Then, generate results for each country:

###
### Afghanistan.
###

# Direct Fire:

# Check the normality of the errors:
panel_af$threshold <- ifelse(panel_af$temperature > 75, 1, 0)

test1 <- lm(df ~ temp_dev*threshold + precipitation + idf + ied + daylength + 
              temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
              precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(DISTID), data = panel_af)
test2 <- lm(log(df+1) ~ temp_dev*threshold + precipitation + idf + ied + daylength + 
              temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
              precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(DISTID), data = panel_af)

test1.stdres = rstandard(test1)
qqnorm(test1.stdres)
qqline(test1.stdres)

test2.stdres = rstandard(test2)
qqnorm(test2.stdres)
qqline(test2.stdres)

temperature_deviation_afghanistan <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temperature > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(df ~ temp_dev*threshold + precipitation + idf + ied + daylength + 
                            temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                            precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                            df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + DISTID | 0 | DISTID, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(df ~ temp_dev*threshold + precipitation + idf + ied + daylength |
                                 week + DISTID | 0 | DISTID, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    
    reg_temp_month <- felm(df ~ temp_dev*threshold + precipitation + idf + ied + daylength +
                             temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                             precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                             df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + DISTID | 0 | DISTID, data = temp)
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(df ~ temp_dev*threshold + precipitation + idf + ied + daylength |
                                  month + DISTID | 0 | DISTID, data = temp)
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    rm(temp, reg_temp_week, reg_temp_week_pars, reg_temp_month, reg_temp_month_pars)
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95 <- temperature_deviation_afghanistan(panel_af, 0.95)
af_90 <- temperature_deviation_afghanistan(panel_af, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
af_95$change <- (sd(panel_af$temp_dev) * af_95$estimate) / mean(panel_af$df)

# For plotting Afghanistan results:
temp_ranges <- unique(af_95$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges)*2-2), by=2)){
  temp_ranges <- append(temp_ranges, " ", after=i)
}

p_af <- ggplot() + 
  
  geom_line(data = af_95[af_95$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = af_95[af_95$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95[af_95$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95[af_95$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95[af_95$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95[af_95$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90[af_90$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90[af_90$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95[af_95$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95[af_95$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90[af_90$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90[af_90$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95[af_95$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95[af_95$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90[af_90$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90[af_90$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95[af_95$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = af_95[af_95$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = af_90[af_90$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = af_90[af_90$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.0002, 0.0008) + 
  annotate(geom = "text", x = 75:105, y = -0.0001825,
           label = temp_ranges,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 45, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# IED:

temperature_deviation_afghanistan_i <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temperature > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + df + idf + daylength + 
                            temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                            precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                            df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + DISTID | 0 | DISTID, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + df + idf + daylength |
                                 week + DISTID | 0 | DISTID, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    
    reg_temp_month <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + df + idf + daylength +
                             temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                             precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                             df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + DISTID | 0 | DISTID, data = temp)
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + df + idf + daylength |
                                  month + DISTID | 0 | DISTID, data = temp)
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    rm(temp, reg_temp_week, reg_temp_week_pars, reg_temp_month, reg_temp_month_pars)
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_i <- temperature_deviation_afghanistan_i(panel_af, 0.95)
af_90_i <- temperature_deviation_afghanistan_i(panel_af, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
af_95_i$change <- (sd(panel_af$temp_dev) * af_95_i$estimate) / mean(panel_af$df)

p_af_i <- ggplot() + 
  
  geom_line(data = af_95_i[af_95_i$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = af_95_i[af_95_i$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_i[af_95_i$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_i[af_95_i$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_i[af_95_i$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_i[af_95_i$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_i[af_90_i$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_i[af_90_i$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_i[af_95_i$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_i[af_95_i$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_i[af_90_i$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_i[af_90_i$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_i[af_95_i$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_i[af_95_i$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_i[af_90_i$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_i[af_90_i$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_i[af_95_i$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = af_95_i[af_95_i$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = af_90_i[af_90_i$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = af_90_i[af_90_i$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.0002, 0.0008) + 
  annotate(geom = "text", x = 75:105, y = -0.0001825,
           label = temp_ranges,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 45, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

afghanistan_wide_combined <- ggarrange(p_af + xlab("") + ylab("OLS"),
                                       p_af_i + xlab("") + ylab(""), 
                                       ncol = 2, nrow=1, labels  = "AUTO")

# PRODUCES FIGURE 9 (PART A):
annotate_figure(afghanistan_wide_combined,
                bottom = text_grob(expression(gamma*" (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### Iraq.
###

# Unconstrained:

temperature_deviation_iraq <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temperature > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(unconstrained1+1) ~ temp_dev*threshold + precipitation + constrained1 + ied + daylength +
                            temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                            precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                            unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                            constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 |
                            week + ADM3NAME | 0 | ADM3NAME, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(unconstrained1+1) ~ temp_dev*threshold + precipitation + constrained1 + ied + daylength |
                                 week + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    
    reg_temp_month <- felm(log(unconstrained1+1) ~ temp_dev*threshold + precipitation + constrained1 + ied + daylength +
                             temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                             precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                             unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                             constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             
                             month + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(unconstrained1+1) ~ temp_dev*threshold + precipitation + constrained1 + ied + daylength |
                                  month + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    rm(temp, reg_temp_week, reg_temp_week_pars, reg_temp_month, reg_temp_month_pars)
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95 <- temperature_deviation_iraq(panel_iq, 0.95)
iq_90 <- temperature_deviation_iraq(panel_iq, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95$change <- (sd(panel_iq$temperature) * iq_95$estimate) / mean(log(panel_iq$unconstrained1+1))
iq_95$change <- (sd(panel_iq$temperature) * iq_95$estimate) / mean(panel_iq$unconstrained1)

# For plotting Iraq results:
temp_ranges <- unique(iq_95$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges)*2-2), by=2)){
  temp_ranges <- append(temp_ranges, " ", after=i)
}

p_iq <- ggplot() + 
  
  geom_line(data = iq_95[iq_95$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95[iq_95$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95[iq_95$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95[iq_95$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95[iq_95$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95[iq_95$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90[iq_90$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90[iq_90$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95[iq_95$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95[iq_95$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90[iq_90$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90[iq_90$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95[iq_95$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95[iq_95$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90[iq_90$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90[iq_90$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95[iq_95$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95[iq_95$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90[iq_90$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90[iq_90$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.0005, 0.0015) + 
  annotate(geom = "text", x = 75:105, y = -0.000475,
           label = temp_ranges,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 45, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# Constrained:

temperature_deviation_iraq_c <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temperature > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(constrained1+1) ~ temp_dev*threshold + precipitation + unconstrained1 + ied + daylength +
                            temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                            precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                            unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                            constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 |
                            week + ADM3NAME | 0 | ADM3NAME, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(constrained1+1) ~ temp_dev*threshold + precipitation + unconstrained1 + ied + daylength |
                                 week + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(constrained1+1) ~ temp_dev*threshold + precipitation + unconstrained1 + ied + daylength +
                             temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                             precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                             unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                             constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 |
                             
                             month + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(constrained1+1) ~ temp_dev*threshold + precipitation + unconstrained1 + ied + daylength |
                                  month + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    rm(temp, reg_temp_week, reg_temp_week_pars, reg_temp_month, reg_temp_month_pars)
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_c <- temperature_deviation_iraq_c(panel_iq, 0.95)
iq_90_c <- temperature_deviation_iraq_c(panel_iq, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_c$change <- (sd(panel_iq$temperature) * iq_95_c$estimate) / mean(log(panel_iq$constrained1+1))
iq_95_c$change <- (sd(panel_iq$temperature) * iq_95_c$estimate) / mean(panel_iq$constrained1)

p_iq_c <- ggplot() + 
  
  geom_line(data = iq_95_c[iq_95_c$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_c[iq_95_c$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_c[iq_95_c$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_c[iq_95_c$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_c[iq_95_c$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_c[iq_95_c$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_c[iq_90_c$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_c[iq_90_c$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_c[iq_95_c$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_c[iq_95_c$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_c[iq_90_c$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_c[iq_90_c$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_c[iq_95_c$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_c[iq_95_c$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_c[iq_90_c$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_c[iq_90_c$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_c[iq_95_c$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_c[iq_95_c$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_c[iq_90_c$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_c[iq_90_c$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.0005, 0.0015) + 
  annotate(geom = "text", x = 75:105, y = -0.000475,
           label = temp_ranges,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 45, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# IED:

temperature_deviation_iraq_i <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temperature > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + unconstrained1 + constrained1 + daylength +
                            temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                            precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                            unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                            constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 |
                            week + ADM3NAME | 0 | ADM3NAME, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + unconstrained1 + constrained1 + daylength |
                                 week + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + unconstrained1 + constrained1 + daylength +
                             temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                             precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                             unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                             constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 |
                             
                             month + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(ied+1) ~ temp_dev*threshold + precipitation + unconstrained1 + constrained1 + daylength |
                                  month + ADM3NAME | 0 | ADM3NAME, data = temp)
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "cluster", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temperature)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temperature)(i + (min(temperature_range)-1))
    
    rm(temp, reg_temp_week, reg_temp_week_pars, reg_temp_month, reg_temp_month_pars)
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_i <- temperature_deviation_iraq_i(panel_iq, 0.95)
iq_90_i <- temperature_deviation_iraq_i(panel_iq, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_i$change <- (sd(panel_iq$temperature) * iq_95_i$estimate) / mean(log(panel_iq$ied+1))
iq_95_i$change <- (sd(panel_iq$temperature) * iq_95_i$estimate) / mean(panel_iq$ied)

p_iq_i <- ggplot() + 
  
  geom_line(data = iq_95_i[iq_95_i$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_i[iq_95_i$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_i[iq_95_i$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_i[iq_95_i$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_i[iq_95_i$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_i[iq_95_i$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_i[iq_90_i$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_i[iq_90_i$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_i[iq_95_i$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_i[iq_95_i$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_i[iq_90_i$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_i[iq_90_i$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_i[iq_95_i$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_i[iq_95_i$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_i[iq_90_i$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_i[iq_90_i$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_i[iq_95_i$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_i[iq_95_i$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_i[iq_90_i$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_i[iq_90_i$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.0005, 0.0015) + 
  annotate(geom = "text", x = 75:105, y = -0.000475,
           label = temp_ranges,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 45, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

iraq_wide_combined <- ggarrange(p_iq + xlab("") + ylab("OLS"),
                                p_iq_c + xlab("") + ylab(""), 
                                p_iq_i + xlab("") + ylab(""), 
                                ncol = 3, nrow=1, labels  = "AUTO")

# PRODUCES FIGURE 9 (PART B):
annotate_figure(iraq_wide_combined,
                bottom = text_grob(expression(gamma*" (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### Create example of process:
###

af_ex <- ggplot() + 
  geom_line(data = panel_af[panel_af$DIST_34_NA %in% "Lashkar Gah",], aes(x = date, y = temperature), color = "gray70") + 
  geom_line(data = panel_af[panel_af$DIST_34_NA %in% "Lashkar Gah",], aes(x = date, y = temp_hat), alpha = 0.5, size = 1.25, color = "blue4") +
  geom_line(data = panel_af[panel_af$DIST_34_NA %in% "Lashkar Gah",], aes(x = date, y = temp_dev), color = "green3") + 
  theme_bw() +
  ggtitle("Lashkar Gah, Afghanistan") + 
  ylab("Temperature (°F)") + 
  xlab("Year") +
  theme(plot.title = element_text(face="bold")) +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold",
                                   size=14),
        axis.title=element_text(size=14,face="bold"))

iq_ex <- ggplot() + 
  geom_line(data = panel_iq[panel_iq$ADM3NAME %in% "Basrah" , ], aes(x = date, y = temperature), color = "gray70") + 
  geom_line(data = panel_iq[panel_iq$ADM3NAME %in% "Basrah" , ], aes(x = date, y = temp_hat), alpha = 0.5, size = 1.25, color = "blue4") +
  geom_line(data = panel_iq[panel_iq$ADM3NAME %in% "Basrah" , ], aes(x = date, y = temp_dev), color = "green3") + 
  theme_bw() + 
  ggtitle("Basrah, Iraq") + 
  ylab("Temperature (°F)") + 
  xlab("Year") +
  theme(plot.title = element_text(face="bold")) +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold",
                                   size=14),
        axis.title=element_text(size=14,face="bold"))

# PRODUCES FIGURE 7:
ggarrange(af_ex, iq_ex)

###
### End.
###
